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ABSTRACT 

A semi-empirical  mathematical  model  for  a turbulent  near  wake  has  been 
developed  within  the  framework  of  an  "interaction  hypothesis”  suggested  by 
Bradshaw'.  The  near  wake  behind  an  airfoil  has  been  treated  as  a "complex" 
shear  flow  consisting  of  two  neighboring  simple  shear  layers  with  dir;tinct  but 
overlapping  shear  stress  profiles  of  opposite  signs.  The  present  model  utilizes 
the  mean  momentum,  and  continuity  equations  together  with  two  shear  stress 
transport  equations  derived  from  the  turbulent  kinetic  energy  equation.  By 
relating  the  shear  stresses  to  the  local  turbulence  quantities,  closure  for 
the  governing  systems  is  achieved  without  the  use  of  the  eddy  viscosity  concept. 
The  shear  stress  is  therefore  no  lonqer  required  to  vanish  at  the  velocity 
extremum.  The  model  has  been  com.parod  to  the  experimental  data  of  Chevray  and 
Kovasznay  and  good  agreement  has  been  obtained.  A comr>arison  with  the  asymmetric 
cascade  wake  data  of  Paj  and  Lakshminarayana  is  also  presented.  Our  calculations 
have  validated  the  basic  philosophy  of  an  interactive  approach  in  the  study  of 
near  wakes.  However,  it  also  clearly  demonstrates  that  the  accuracy  of  certain 
empirical  functions  used  to  define  the  turbulence  structure  has  a direct  impact 
on  the  success  of  this  or  any  other  calculation  method. 
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INTPODL'CTTON 

One  of  the  primary  objectives  of  the  present  research  program  is  to  develop 
an  accurate  and  effective  technique  for  predicting  aerodynamical ly  induced  vibra- 
tions of  rotor  and  stator  blades  in  turbomachines.  From  a theoretical  standpoint, 
an  essential  step  towards  the  development  of  such  a method  requires  a realistic 
modelling  of  the  flow  through  an  annular  cascade  of  airfoils.  resulting  prob- 

lem is,  of  course,  three-dimensional,  but  it  can  be  substantially  simplified  if 
we  assume  axial  symmetry  in  the  time-averaaed  basic  flow  so  that  a two-dimc.-nsional 
analysis  in  the  circumferential  or  blade-to-blade  stream  surface  is  possible. 

Within  this  general  context,  the  turbulent  near  wake  behind  any  one  airfoil  con- 
stitutes but  a small  part  of  the  overall  cascade  flow  field.  Nevertheless,  the 
structure  of  a near  wake  plays  a decisive  role  in  determ.ining  the  aerodynamic 
forcing  function  which  will  be  needed  for  any  subsequent  cascade  gust  calculations. 

Despite  recent  advances  in  numerical  methods  dealing  with  turbulent 
boundary  layers,^  little  progress  has  been  made  in  the  modelling  of  an  asymmetric- 
turbulent  near  wake.  This  is  in  part  due  to  the  fact  that  the  state  of  the  art 
in  turbulent  modelling  is  such  that  a large  amount  of  empiricism  is  often  required 
to  insure  the  success  of  any  prediction  method.  In  the  case  of  a near  wake  formed 
by  the  coalescence  of  two  turbulent  boundary  layers  with  opposite  shear  stresses, 
very  little  is  known  cibout  the  mechanism  of  their  interaction.  Moreover,  for  an 
asymmetric  wake,  the  point  of  zero  Reynolds  stress  does  not  necessarily  coincide 
with  the  point  of  zero  velocity  gradient.  Any  attempt  to  obtain  closure  for  the 
governing  equations  by  using  a simple  eddy  viscosity  model  for  the  shear  stress 
will  therefore  encounter  difficulties.  It  is  the  purpose  of  this  work  then  to  out- 
line a method  by  which  these  problems  can  be  resolved. 

In  important  element  in  our  present  analysis  is  the  "interaction  hypothesis" 

2 

first  proposed  by  Bradshaw,  et.  al.  in  conjunction  with  the  study  of  duct  flow. 
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within  the  context  of  this  hypothesis,  a near  wake  can  be  reaarded  as  an  inter- 
action between  two  neighboring  simple  shear  layers  with  distinct  but  possibly 

i 

overlapping  shear  stress  profiles  of  opposite  signs  (cf.  Figure  1).  If  the  j 

interaction  is  such  that  the  tuAbuZ^nce  ■itAuctuAe.  in  each  layer  is  essentially 
unaffected  by  the  presence  of  the  adjacent  layer,  a superposition  of  the  two  shear 
stress  fields  for  the  purpose  of  calculating  the  net  Reynolds  stress  is  then  | 

possible.  Moreover,  by  relating  the  shear  stresses  to  the  local  turbulent  quanti-  | 

ties,  as  it  has  been  suggested  by  Bradshaw,  Ferriss  and  Atwell,  two  shear  stress  i 

transport  equations  can  be  derived  from  the  turbulent  kinetic  energy  equation. 

Closure  for  the  governing  system  is  therefore  achieved  without  making  use  of  the 
eddy  viscosity  concept. 


I 


FIGURE  1.  THE  INTERACTION  HYPOTHESIS 


BASIC  FOR-VVIJiTION:  INTERACTIVE  APPROACH 


In  order  to  explore  the  validity  of  the  interaction  approach,  our  inver.ti- 
gation  has  been  confined  to  a prototype  model  which  includes  only  the  essential 
physical  details.  The  fluid  is  initially  assumed  to  be  incompressible.  Further- 
more, due  to  the  magnitude  of  the  Reynolds  number  usually  encountered  in  turbo- 
machine/internal flow  applications,  the  flow  is  also  expected  to  satisfy  the 
boundarv’  layer  appro.ximation  and  t.he  viscous  stresses  have  been  ignored  when  com- 
pared to  the  Reynolds  stresses.  The  momentum  and  continuity  equations  are  then 
given  by 
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The  Cartesian  coordinate  system  for  those  equations  is  chosen  such  that  the 
x-axis  is  tangent  to  the  mean  camb( r line  of  an  upstream  airfoil  at  the  trailing 
edge;  x is  equal  to  0 at  the  trailing  edge  and  positive  in  the  downstrear  direc- 
tion. The  time-averaged  velocity  com.ponents  along  the  x and  y axes  are  then  given 
by  U and  V,  respectively.  The  free  stream  velocity  in  the  pressure  gradient 
term  of  equation  (1)  is  to  be  replaced  by  L'^^or  U depending  on  whether  y > y^ 

or  y < y , where  y (x)  is  the  locus  of  the  velocity  minimum.  Physically,  our 
c c 

prototype  problem  can  be  envisioned  as  that  of  a two-dimensional  asymmetric  wake 
behind  an  airfoil  or  a flat  plate,  formed  by  the  coalescence  of  the  upper  and 
lower  surface  boundary  layers  having  possibly  different  characteristics.  In  our 
present  notation,  the  usual  shear  stress  is  replaced  by  T = (shear  stress)/ 
(density)  = - uv,  where  u and  v are  velocity  fluctuations  and  the  overbar  denotes 


a time  average. 


r. 


According  to  the  interaction  hypothesis,  wo  can  further  write 

+ 

T = T + T 

such  that  is  the  dominant  shear  stress  in  the  region  of  positive  shear — 

nominally  where  y > y^ — and  vice  versa  for  x . In  order  to  obtain  a tractable 

system,  certain  assumptions  must  be  made  to  secure  closure  for  equations  (1) 

3 

and  (2).  Using  the  Bradshaw-Ferriss-Atwell  turbulence  model  we  define  two 

^ — 

sets  of  empirical  functions,  (a^  ,G  ) and  (a^  ,G  ,L  ),  such  that 

+ 

± T 
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(3) 
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(p'v  + 
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q^v) 

l/2~ 


(4) 


and 


1/2 


L"  = . (5) 

The  relations  (3)  - (51  are  equivalent  to  assuming  that  (i)  the  local  shear 
stress  is  proportional  to  the  turbulent  intensity  ; (ii)  the  energy  diffu- 
sion is  directly  proportional  to  the  local  shear  stress  with  a factor  depending 

on  the  maximum  of  th^-  shear  stress,  x , and  (iii)  the  dissipation  rate  C is 

m 

determined  by  local  shear  stress  and  a length  scale  L.  With  equations  (3),  (4) 

and  (5),  the  structure  of  turbulence  is  defined  in  terms  of  relations  between 

turbulence  quantities  and  it  is  thus  independent  of  the  mean  flow.  It  has  been 

observed  that  the  above  assumptions  are  quite  valid  over  a wide  range  of  pressure 

gradients  V.ach  number,  Reynolds  number  and  appear  to  be  rather  insensitive  to 

4 , 5 

rapid  chanoe.s  in  the  mean  flow.  This  will,  of  course,  enable  the  present 

method  to  be  extended  to  a wide  variety  of  mean  flow  conditions,  often  without 
the  need  tor  extra  data.  On  substituting  (3),  (4)  and  (5)  into  the  exact  turbu- 
lent energy  equation,  we  can  derive  two  transport  equations  for  X^  and  x of  the 
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+ + 1/2 


/„  9 . „ 3 V X \ _±  T 3 /^-X  ±1  i|l/2\  - 1 i 

^ 9^)(  — T ’ = " " 37  r ^ '"-X  i " ^ 


The  system  consistinq  of  equations  (1) , (2)  and  (6)  is  now  closed  provided 

that  the  empirical  functions  can  be  determined  experimentally.  A schematic 
representation  of  these  functions  is  qiven  in  Figure  2.  Conceptually,  the  two 
sets  of  empirical  functions  are  to  be  regarded  as  distinct,  but  in  fact  they 
are  qualitatively  the  same  except  for  changes  in  signs  and  their  dependence  on 
y as  illustrated  in  Figures  2 and  3.  It  must  be  emphasized  that  the  success 
of  the  present  approach,  regardless  of  the  method  used  in  solving  the  governing 


system  of  equations  is  then  dependent  on  the  accuracy  of  (a^~,G' ,L~).  The 
details  of  the  empirical  functions  used  in  our  calculations  are  discussed  in 
the  RESULTS  AND  DISCUSSION  section.  In  addition  to  equations  (1),  (2)  and  (6), 


certain  initial  and  boundary  conditions  m.ust  be  specified. 

On  noting  that  the  genesis  of  a near  wake  is  effected  by  the  merging  of 
two  surface  boundary  layers  at  the  trailinq  edge  where  x = 0,  the  natural  initial 
conditions  for  an  airfoil  with  zero  trailinq  edge  thickness  are  therefore  qiven 


U = U(0,y)  , 


V - V(0,y)  , 


^ i 0 


y > 0 I 
y < 0 i 


, y > 0 


T (0,y)  , y < 0 


UPPER  WAK 


FIGURE  2.  A SCHEMATIC  REPRESENTATION  OF  THE  EMPIRICAL  FUNCTIONS 
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where  the  velocity  profiles  U(C,y)  and  V(0,y),  and  the  shear  stress  profiles 
T (0,y)  and  T (0,y)  are  ideally  supplied  by  an  upstream  boundary  layer  calcu- 
lation or  measured  directly  from  experiments.  For  x > 0,  we  shall  also  require 
that  V(x,y^)  = 0,  when  y^  = 0 for  the  symmetric  case. 

The  boundary  conditions  for  the  present  problem  are 


u -*■  u 


as  y -*•  + °°  , 


(11) 


and 


as  y 


(12) 


A careful  examination  of  the  turbulent  shear  stress  equations  also  indicates 
that  should  vanish  almost  immediately  after  crossino  into  the  reaion  domi- 
nated by  negative  shear.  In  fact,  on  noting  that  the  cjoverning  system  is  hyper- 
bolic,  it  is  p'ossible  to  determine  two  characteristics,  . and  , , emanating 

from,  the  trailing  edge  such  that  t"*”  = 0 in  a region  bounded  by  T and  the  lower 
edge  of  the  wake,  and  vice  versa  for  T . In  differential  form,  these  characteris- 
tics are  given  by 

r~ : ^ = tan  V and  y(0)  = 0 , (n) 

dx  ± 


where 


tan  = 


V -t  a/  G“  T 
1 ' m 


*+,  + + i . . 

t [(a.,  n-)’!T  + 2a,  1 ] 

1 m 1 


U 


Hyperbolicity  has  indeed  been  fully  exploited  in  the  work  of  Pradshaw,  et . 
al.(see,  e.g.,  references  2,  3 and  4).  In  the  present  analysis,  however,  the 
characteristics  associated  with  the  governing  system  play  a relatively  minor 
role.  They  are  useful  in  so  far  as  providing  a conceptual  framework  within 
which  the  precise  boundaries  of  the  interacting  shear  layers  can  be  defined  for 


■f 


the  purpose  of  scalinq  the  empirical  functions  (cf.  Fiqures  2 and  3) . On  the 
other  hand,  the  conditions  that  and  T be  respectively  zero  in  most  of  th*-’ 
reqions  of  neqative  and  positive  shear  will  be  automatically  satisfied  if  we 
require  that 


T-*Oasy-*-_oo  (14) 


and  T-+Oasy-*+*.  (15) 

The  explicit  use  of  characteristics  in  specifyina  boundary  conditions  is  there 
fore  unnecessary.  It  is  useful  to  note  here  a few  key  points  of  the  interacti 
analysis : 


(i)  Due  to  the  inclusion  of  the  boundary  layer  profiles 
at  the  trailino  edae  as  initial  conditions,  the 
observed  flow  asymmetry  is  beinq  dealt  with  directly. 

(ii)  Since  the  near  wake  is  hiahly  dependent  upon  the 
structure  of  the  lower  and  upper  surface  boundary 
layers,  the  interaction  approach  is  expected  to 
produce  better  results  than  any  similarity  methods 
in  which  the  loss  of  initial  conditions  is  implicit 
in  the  zero  initial  wake-width  assumption. 

(iii)  The  (indirect)  use  of  the  turbulent  enerqy  equation  loads 
to  the  incorporation  of  advective  and  diffusive  effects 
in  the  enerqy  budqet , in  addition  to  the  usual  balance 
between  production  and  dissipation  of  turbulence.  The 
"past  history"  of  the  turbulence  is  therefore  explicitly 
taken  into  account.  Tn  a hiqhly  turbulent  and  rapidly 
varv'inq  flow  reqion,  such  as  a near  wake,  this  considera- 
tion will  have  a direct  bearinq  on  the  accuracy  of  the 


present  model. 
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(iv)  The  decomposition  of  the  net  shear  stress  profile  into 
+ 

and  T , each  of  which  satisfies  a separate  transfxjrt 
equation,  enables  us  to  circumvent  the  usual  difficulties 
associated  with  any  complex  shear  tiow  havina  an  extremum 
in  its  velocity  profile. 

(v)  It  is  also  of  interest  to  point  out  that  past  experience 
with  boundary  layers  and  duct  flow  indicates  that  the 
empirical  functions  exhibit  a remarkable  ne-asure  of 
universality.^  It  is  exi  ected  that  the  corresponding 
functions  for  near  wakes,  once  determined,  will  be  canon- 
ical, to  some  extent,  to  the  present  c! ass  of  flow.  In 
any  internal  flow  applications,  a particular  cascade 
geometry  and  the  angle  of  attack  or  incidence  of  the 
airfoil  will  only  influence  the  near  wake  solution 
through  the  initial  and  boundary  conditions. 


In  the  foreaoing  analysis,  the  fluid  has  been  assumed  to  be  incompressible. 
This  simplification  can  be  removed  without  added  complexities  if  we  invoke 
Morkovin's  hy{X)thesis  that  the  turbulence  structure  — defined  in  our  case 
by  the  empirical  functions  --  is  unaffected  by  compressibility  when  the  Mach 
number  fluctuation  is  much  less  than  unity.  This  approach  has  been  found 
to  be  quite  satisfactory  by  Bradshaw^  for  boundary  layers  on  adiabatic 
walls  at  free-stream  Mach  numbers  up  to  4.  Thus,  using  an  order  of  magnitude 
analysis  identical  to  that  of  Bradshaw,  the  mean  momentum,  continuity,  and 
turbulent  kinetic  energy  equations  for  adiabatic  compressible  flow  take  on 
the  following  simplified  forms: 


21 

9x 


(V  + 


p' V 9u  _ _ J:_  ^ J_  5pi 

0 9y  — 9x  — 9y 

D 0 


(16) 


14 


9U 


* r(Y 


4- 


?U  . V 
3y  U 


r(Y  - DM  ^ 
00 


U dp  _ 
— dx 


P 


(25) 


U 


DM  ^ 

OO 


•f  4- 

T__  9t 

U 3y 


■ 

4- 

1 


r (Y 


1) 


iH 


+ 


?y 


(G"t 


(26)  - (27) 


where  the  coefficient  in  the  pressure  qradient  term  of  equation  (24)  has 
been  rewritten  such  that  the  mean  density  does  not  appear  ex^'licitly.  The 
free  stream  Mach  number  M^  in  equations  (24)  and  (25)  must  be  replaced  by 
M^  or  M depending  on  the  shear  region,  and  the  pressure  gradient  must 
also  be  adjusted  accordinoly . It  can  easily  be  shown  that  equations  (24j  - 
(27)  are  again  hyperbolic  and  they  are  reducible  to  the  incompressible 
equations  when  M 0 . The  solution  algorithm  to  be  discussed  in  the 
next  section  is  therefore  equally  applicable  to  the  compressible  case 
with  only  minor  modifications. 
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SOLUTION  ALGORITHM 


A n'jinber  of  very  effective  alqorithms  have  been  developed  in  recent  years 
for  solving  quasi-linear  hyperbolic  systems  of  the  form  given  by  equations 
(1) , (2)  and  (6).  In  the  work  of  Bradshaw,  et.al.  on  turbulent  boundary  layers, 
the  method  of  characteristics  was  adopted  to  solve  a somewhat  simpler  system,  of 
three  equations  involvino  the  dependent  variables  U,  V and  T.  The  advantages 
of  such  an  apcroach  are  that  the  resulting  numerical  scheme  is  essentially 

explicit,  and  that  the  well-known  Courant-Friedr ichs-Lewy  criterion  provides  a 
clear  cut  limit  to  the  x-step  that  can  be  taken  to  insure  numerical  stability. 

The  uncoupling  of  the  continuity  equation  from  the  remaining  system  in  that  case 
also  effectively  reduces  the  problem  to  one  of  only  two  unknowns — namely,  U and 
r--so  that  the  characteristic  angles  at  each  point  in  the  flow  can  be  con- 
veniently determined  by  a second-order  algebraic  equation. 

For  our  present  purpose,  a direct  generalization  of  the  method  of 
characteristics  to  deal  with  equations  (1),  (2)  and  (6)  is,  of  course,  pos- 
sible but  has  been  found  to  be  somewhat  cumbersome.  Despite  the  fact  that 
V can  again  be  uncoupled,  the  determination  of  the  characteristic  angles  now 
requires  the  solution  of  essentially  a third-order  algebraic  equation.  An 

attempt  to  circumvent  this  difficulty  has  since  been  discussed  by  Bradshaw, 

2 

Dean  and  McEligot  in  conjunction  with  their  work  on  duct  flow.*"  Due  to  the 
rather  weak  coupling  of  the  governing  system,  it  has  been  suagested  that  the 
two  turbulent  shear  stress  equations  can  be  solved  separately  using  the 
velocity  profile  at  the  previous  x-step  as  a first  approximation.  The  result- 
inq  T and  T are  then  summed  and  a new  velocity  profile  is  then  calculated. 
Accuracy  can  presumably  be  improved  throuQh  iteration.  For  our  purpose,  however 
we  find  it  convenient  to  formulate  a more  direct  numerical  scheme  based  on  the 

7 

general  finite  difference  approach  similar  to  the  one  devised  by  Ferriss. 


If) 


In  anticipation  of  the  rapidly  varyina  velocity  profiles  in  the  inner 
wake  reaion  similar  to  those  encountered  in  boundary  layers,  a mesh  of  variable 
qrid  si?.e  in  the  normal  directions  has  been  adopted  for  the  present  analysis. 
The  y-steps  are  qradated  outward  from  the  x-axis  in  the  form  of  a geometric 
sequence  such  that 

I j ! 

(Ay)^  = a(Ay)^_^  = a (Ay)|^, 

where 

(28) 

= y.,,-  y,  . 

and 

(Ay)o 


The  geometric  ratio  u is  usually  taken  to  be  of  the  order  of  1.05.  We  shall 

also  denote  the  boundary  qrid  points  for  the  lower  and  upper  wake  edaes 

respectively  by  y and  v,  . The  total  number  of  interior  grid  points  at  any 
-L  M 

given  x-station  is  then  given  by 

J = L + M - 1 . 


The  required  numerical  scheme  for  solvina  equations  (1),  (2)  and  (6)  can 
now  be  formulated  in  terms  of  difference  quotients  witnin  a typical  computa- 
tion cell  consisting  of  two  adjacent  mesh  rectangles  as  shown  in  Figure  4. 

For  any  function  f(x,y)  which  is  at  least  twice  differentiable  in  /,  the  second- 
order-correct  finite  difference  approximation  to  the  partial  derivative  of  f 
with  respect  to  y at  (x^,y^)  is  then  given  by 


where 


• (l->^) 
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ill) 
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j < 1 I 

(29) 
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In  the  case  of  a logarithmic  function  of  y,  it  can  be  shown  that  the  error 
associated  with  an  approximation  of  the  form  given  by  (29)  remains  constant 
throughout  the  entire  grid  despite  the  increase  in  y-step  size  in  the  outer 
regions  (e.g.,  reference  8).  A"centered" difference  equation  is  used  to  approxi- 
mate the  partial  derivative  of  f with  respect  to  x.  Thus  we  have 


irJ. 


i+1/2,  j 


f,  , f.  , 

1,1 

(Ax)  . 

1 


where 


’'i  + 1/2 


X.  , + X. 

1 + 1 1 


(30) 


and 


(Ax) . = X.  , - X. 
1 1 + 1 1 


Moreover,  it  is  useful  to  note  that 
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The  finite  difference  analog  to  the  governing  system  can  now  be  derived 
by  replacing  the  partial  derivatives  in  equations  (1),  (2)  and  (6)  with  their 
corresponding  approximations  given  by  (30)  and  (31) . For  this  pur[X3se,  we  now 
de fine 
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where 


+ 

e7 


V.  , ^ , 1/2 

= ± G-  It-  ! 

i,D,r  ± i+l,j+r'  m'.  r = -1,  0,  1, 

2ai  1 


and 


(&Y) 


.,j  (Ax)^ 


- (A+1)  u.  . 

1.1 


Tt  can  then  be  shown  that  around  the  point  j ^ 

and  (6)  are  reducible  to 
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The  inhomogeneous  term  d^  ^ is  given  by 
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In  equation  (42),  the  edqe  velocity  U must  be  replaced  by  U or  U , 

‘CO  ^ ‘ 4.00  — on' 

depending  on  whether  y^  > y^  or  y^  ^ course  of  deriving  equation 

( 36)  , we  have  also  linearized  the  corresponding  equations  by  approximating  the 
coefficients  involving  U,  V,  or  T by  their  values  at  x^,  which  are  ta)cen  to 
be  known.  The  non-linear  version  m.ay  be  recovered,  if  needed,  by  replacing 
i + 1/2  for  the  i-indices  of  those  dependent  variables  appearing  in  equations 
(32)  - (35).  At  any  given  x-statiop,  the  J interior  arid  points  then  generate 
3J  equations  with  the  same  number  of  unltnowns.  The  resulting  system  can  be 


written  in  the  following  compact  form: 
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where  the  i- indices  have  been  omitted  for  clarity.  The  boundary  conditions  at 

y _ and  y„  are  incorporated  throuqh  the  inhomoqeneous  vector  d.  Before  the  above 
-L  M 

linear  system  is  actually  solved,  however,  a modification  is  desirable.  On  notinq 
that  the  downstream  shear  stress  distributions  of  both  the  top  and  bottom  shear 
layers  should  vanish  in  most  of  the  regions  of  opp'ositc  shear  due  to  zero  up- 

4-  — 

stream  conditions,  it  is  possible  to  calculate  T and  T only  at  those  grid 
points  where  they  are  non-zero.  The  advantage  of  this  approach  is  that  it  re- 
sults in  a reduction  of  the  number  of  equations  and  therefore  the  size  of  the 
matrix  M.  A more  important  consequence  is  that  we  are  dispensed  with  the  need 
to  define  and  justify  the  empirical  functions  too  far  beyond  the  boundaries 
within  which  the  top  and  lower  shear  layers  are  confined.  On  the  other  hand, 
care  must  be  exercised  to  allow  for  the  growth  of  the  shear  layers.  For  tViis 
purpose,  whenever  a shear  stress  profile  is  computed,  the  magnitude  of  the 
shear  stress  at  an  edge  grid  point  is  checked  to  determine  whether  certain  limits 


(usually  0.01  l have  been  exceeded.  The  boundary  for  the  downstream  profile 


will  then  be  moved  outward  by  one  grid  point  if  the  tolerance  has  not  been  met. 

The  reduced  form,  of  the  coefficient  matrix  is  of  the  band  type  with 
nine  non-zero  diagonals.  The  liiiear  system  can  therefore  be  solved  efficiently 
by  one  of  the  many  standard  routines  available  through  various  scientific  sub- 

4-  _ 

routine  packages.  The  values  of  U,  T and  T at  the  interior  grid  points  are 
then  contained  in  the  solution  vector 

Our  linearization  of  the  numerical  scheme  also  ’..as  the  effect  of  un- 
coupling the  continuity  equation,  and  therefore  V,  from  the  remaining  system. 
Using  equation  (2)  to  eliminate  (f^U/J^x)  from  the  momentum  equation,  we  obtain 

dU 


r,2  3 ,V.  „ ""on 

- V 75 — ( — ) - U — + TC- 

3y  u dx  3y 


(44) 


where  all  quantities,  except  V,  are  now  known.  L'sinq  the  initial  condition 

V(y  ) = 0,  equation  (32)  can  then  be  inteqrated  outwards  from  y to  determine  the 
c c 

vertical  velocity  profile  V. 

With  the  values  of  U,  T^,  X and  V known  at  x^,  the  implicit  finite  differ- 
ence scheme  can  be  used  to  advance  U,  x^  and  X to  x . V is  then  calculated 

1 + 1 

by  equation  (44).  This  process  is  repeated  until  a solution  over  the  entire  do- 
main of  interest  is  obtained.  Alternately,  before  an  increment  in  the  x-direction 
is  made,  the  current  velocity  and  shear  stress  profiles  can  be  used  to  update  the 
linearized  coefficients  and  to  improve  the  solution  by  iteration. 

The  proposed  numerical  scheme  is  believed  to  be  quite  stable,  but  caution 

must  be  exercised  to  insure  the  accuracy  of  the  solution.  Although  the  effects 

of  a Goldstein-type  singularity  at  the  trailing  edge  is  exj.>ected  to  vanish 
g 10 

quickly,'  a recent  study  by  Burggraf  suggests  that  an  extremely  fine  mesh 

must  be  used  to  insure  accuracy  in  the  very  near  wake  region.  XYpically  an  x- 

-3/5 

step  size  of  the  order  of  R , R beino  the  Reynolds  numl:>er  based  on  chord 
length,  must  be  used  for  a mixing  length  model.  Our  numerical  experiments 
indicate  a comparable  x-step  size  must  also  be  used  (for  the  first  2 or  3% 
chord  length  downstream)  in  the  present  calculation.  Moreov'er,  it  is  useful 
to  note  that  the  rates  of  growth  for  the  lower  and  upper  wake  edges  can  be 


determined  by  the  maximum  characteristic  anales  y and  Y at  y and  y 
' +0-  -L  M 


From  equation  (13)  we  derive 


tan  Y, 


4 +. 


^ 2a,”G“  X' 

1 ' m' 


Thus,  the  x-step  size  must  also  in  aeneral  be  restricted  such  that  the  wake 
width  is  allowed  to  expand  by  no  more  than  one  arid  point  in  either  normal 
direction  for  each  x-increment. 
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RESULTS  AND  DISCUSSION 

As  the  first  test  case  for  our  rnodel , we  have  chosen  the  measurements  by 
Chevray  and  Kovasznay  of  a symmetric  wake  behind  a thin  flat  plate. The  data 
available  in  this  case  is  perhaps  the  most  extensive  and  best  documented  from 

9 

both  the  experimental  and  the  theoretical  points  of  view.  The  study  by  Burqgraf 

adapted  the  Cebeci-Smith  and  Glushko  models  to  compute  the  velocity  and  shear 

stress  profiles  and  it  thus  provide.s  a basis  of  comparison  between  the  present 

model  and  other  prediction  methods  now  available.  It  should  be  noted,  however, 

that  due  to  the  use  of  the  eddy  viscosity  concept,  the  Cebeci-Smith  and  Glushko 

models  used  by  Burqqraf  are  restricted  to  symmetric  fiows.  The  present  approach 

does  not  assume  symmetry  in  its  formulation  and  a comparison  with  the  asymmetric 

12 

cascade  wake  data  of  Raj  and  Lakshminarayana  is  presented  as  the  second  test 
case . 

As  repeatedly  noted  in  our  previous  discussions,  one  of  the  most  crucial 
steps  in  the  application  of  the  present  scheme--or,  in  fact,  any  other  calcula- 
tion method — is  the  choice  of  empirical  functions  which  define  the  turbulence 
structure.  In  a recent  interactive  study  of  symmetric  jets  and  wakes  by  Morel, 
a set  of  empirical  functions  (see  Figure  5)  were  suggested  for  wake  calculations 
starting  at  a point  well  downstream  of  the  trailing  edge.  The  present  solution 
algorithm  has  been  used  to  reproduce  Morel's  results  as  a check  and  good  agree- 
ment was  obtained;  but  Morel's  functions  are  not  valid  in  the  near  wake  region. 

It  can  be  argued  that  the  evolution  of  two  coalescing  boundary  layers  into 
a near  wake  is  initially  confined  to  a region  bounded  by  the  two  characteristics 
emanating  from  the  trailing  edge.^^  As  a first  approximation , the  boundary  layer 
length  scale  L is  expected  to  be  valid  except  near  the  wake  center-line  where  it 
has  been  assumed  to  be  constant  and  proportional  to  the  width  of  the  "inner  retnon" 


(From  Morel,  1972) 
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(see  Figure  6).  We  have  tentatively  used  Morel's  function  G and  the  ixjint  where 
G = 0 is  correlated  with  the  shear  stress  maximum.  With  this  normalization,  the 
function  G closely  approximates  its  boundary  layer  counterpart  in  the  outer  wake 
edges  but  assumes  negative  values  near  the  wake  center-line.  This  allows  the 
shear  stresses  to  be  diffused  away  from  their  maxima.  The  value  for  a^  is  again 
taken  to  be  0.15.  The  empirical  functions  used  in  our  present  calculations  are 
therefore  necessarily  crude  pendina  more  thorouch  experimental  investigations. 

In  particular,  the  effects  of  the  ass’jmed  "interaction"  between  the  two  shear 
layers  must  be  clarified  and  incorporated  as  adjustments  in  the  em,pirical  func- 
tions. Nevertheless,  the  results  obtained  for  the  symmetric  case  are  encouraging 
and  they  are  presented  in  Figures  7-12  along  with  results  obtained  by  Burggraf 
based  on  the  Cebeci-Fmith  and  Glushko  models.  The  present  method  produces  good 
aqreements  with  the  experimental  velocity  profiles  at  = 8.6,  34.4  and  86.2 

--corresporuincr  to  x/c  = 0.0208,  0.083  and  0.208  respectively.  The  predicted 
shear  .stress  ma.xima  also  show  improvements  over  those  of  Cebeci-.Smith  and  Glushko. 

The  present  m.odel  has  also  been  used  to  compute  the  asymmetric  wake  pro- 
files behind  a cascade  of  airfoils  at  an  ancle  of  incidence  i = -6°,  and  com- 
parisons with  the  exnerimental  measurements  of  Raj  and  I.akshminarayana  at  x/c  = 
0.08,  0.24  and  0.40  are  given  in  Figures  13  - 18.  The  empirical  functions  are 
exactly  those  used  in  the  symmetrical  case.  The  inherent  asymmetry  is  accounted 
for  through  the  previously  discussed  scalinc  procedure.  The  overall  acreements 
between  the  predicted  velocity  profiles  and  the  data  is  cood  and  the  present  re- 
sults show  considerable  improvements  over  those  obtained  by  either  the  global  or 

14 ,15 

the  local  similarity  methods  --particularly  at  the  near  wake  stations.  It 

should  be  pointed  out,  however,  that  the  results  for  the  asymmetric  case,  though 
satisfactory , do  not  agree  as  closely  with  the  experimental  results  as  those  for 
the  symmetric  case.  This  is  at  least  in  part  due  to  the  fact  that  the  observed 
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Figure  6.  The  empirical  function 
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Mean  velocity  distribution  in  wake 
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Mean  velocity  distribution  in  wake. 
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Figure  10 

Distribution  of  Raynolds  strass. 
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Figure  1 2 

Distribution  of  Reynolds  stress. 
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turbulent  shear  stress  level  for  the  asyiranetric  cascade  wake  is  much  higher  than 
that  for  the  symmetric  case  and  the  accuracy  of  the  empirical  functions  is  more 
critical.  It  is  also  of  interest  to  note  that  the  accuracy  of  the  present  model 
seems  to  deteriorate  as  one  progresses  downstream  of  the  airfoil  trailing  edge 
(cf.  Figures  16  through  18) . The  simple  adjustments  made  in  the  length  scale 
function  L in  our  attempt  to  account  for  the  evolution  of  two  coalescing  boundary 
layers  into  a fully  developed  far  wake  are  perhaps  too  crude  to  adequately  model 
this  extremely  complex  and  intricate  process. 

In  this  connection,  it  is  highly  desirable  that  some  detailed  experimental 
measurements  with  sufficient  spatial  resolution  be  carried  out  in  the  near  future 
in  order  to  obtain  a good  data  base  from  which  the  boundary  layer/wake  mixing  can 
be  quantified  and  the  physical  basis  for  the  interaction  hypothesis  be  ascertained. 
These  experiments  should  also  result  in  a concomittant  generation  of  accurate 
empirical  functions  which  are  correlated  to  the  basic  properties  of  the  shear  flows. 

In  summary,  our  investigation  so  far  has  validated  the  basic  philosophy  of 
the  interactive  approach  for  near  wake  calculations,  but  it  also  clearly  suggests 
that  further  refinements  in  the  empirical  functions  based  on  reliable  experimental 
data  are  needed  before  the  present  method  can  be  used  with  confidence  as  a design 
tool  for  an  arbitrary,  complex,  turbulent  shear  layer. 


r 
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